Reaction-Based Enumeration, Active Learning, and Free Energy Calculations To Rapidly Explore Synthetically Tractable Chemical Space and Optimize Potency of Cyclin-Dependent Kinase 2 Inhibitors.

The hit-to-lead and lead optimization processes usually involve the design, synthesis, and profiling of thousands of analogs prior to clinical candidate nomination. A hit finding campaign may begin with a virtual screen that explores millions of compounds, if not more. However, this scale of computational profiling is not frequently performed in the hit-to-lead or lead optimization phases of drug discovery. This is likely due to the lack of appropriate computational tools to generate synthetically tractable lead-like compounds in silico, and a lack of computational methods to accurately profile compounds prospectively on a large scale. Recent advances in computational power and methods provide the ability to profile much larger libraries of ligands than previously possible. Herein, we report a new computational technique, referred to as "PathFinder", that uses retrosynthetic analysis followed by combinatorial synthesis to generate novel compounds in synthetically accessible chemical space. In this work, the integration of PathFinder-driven compound generation, cloud-based FEP simulations, and active learning are used to rapidly optimize R-groups, and generate new cores for inhibitors of cyclin-dependent kinase 2 (CDK2). Using this approach, we explored >300 000 ideas, performed >5000 FEP simulations, and identified >100 ligands with a predicted IC50 < 100 nM, including four unique cores. To our knowledge, this is the largest set of FEP calculations disclosed in the literature to date. The rapid turnaround time, and scale of chemical exploration, suggests that this is a useful approach to accelerate the discovery of novel chemical matter in drug discovery campaigns.

WScore: A Flexible and Accurate Treatment of Explicit Water Molecules in Ligand-Receptor Docking.

We have developed a new methodology for protein-ligand docking and scoring, WScore, incorporating a flexible description of explicit water molecules. The locations and thermodynamics of the waters are derived from a WaterMap molecular dynamics simulation. The water structure is employed to provide an atomic level description of ligand and protein desolvation. WScore also contains a detailed model for localized ligand and protein strain energy and integrates an MM-GBSA scoring component with these terms to assess delocalized strain of the complex. Ensemble docking is used to take into account induced fit effects on the receptor conformation, and protein reorganization free energies are assigned via fitting to experimental data. The performance of the method is evaluated for pose prediction, rank ordering of self-docked complexes, and enrichment in virtual screening, using a large data set of PDB complexes and compared with the Glide SP and Glide XP models; significant improvements are obtained.

Improved docking of polypeptides with Glide.

Predicting the binding mode of flexible polypeptides to proteins is an important task that falls outside the domain of applicability of most small molecule and protein-protein docking tools. Here, we test the small molecule flexible ligand docking program Glide on a set of 19 non-α-helical peptides and systematically improve pose prediction accuracy by enhancing Glide sampling for flexible polypeptides. In addition, scoring of the poses was improved by post-processing with physics-based implicit solvent MM-GBSA calculations. Using the best RMSD among the top 10 scoring poses as a metric, the success rate (RMSD ≤ 2.0 Å for the interface backbone atoms) increased from 21% with default Glide SP settings to 58% with the enhanced peptide sampling and scoring protocol in the case of redocking to the native protein structure. This approaches the accuracy of the recently developed Rosetta FlexPepDock method (63% success for these 19 peptides) while being over 100 times faster. Cross-docking was performed for a subset of cases where an unbound receptor structure was available, and in that case, 40% of peptides were docked successfully. We analyze the results and find that the optimized polypeptide protocol is most accurate for extended peptides of limited size and number of formal charges, defining a domain of applicability for this approach.

Docking performance of the glide program as evaluated on the Astex and DUD datasets: a complete set of glide SP results and selected results for a new scoring function integrating WaterMap and glide.

Glide SP mode enrichment results for two preparations of the DUD dataset and native ligand docking RMSDs for two preparations of the Astex dataset are presented. Following a best-practices preparation scheme, an average RMSD of 1.140 Å for native ligand docking with Glide SP is computed. Following the same best-practices preparation scheme for the DUD dataset an average area under the ROC curve (AUC) of 0.80 and average early enrichment via the ROC (0.1 %) metric of 0.12 were observed. 74 and 56 % of the 39 best-practices prepared targets showed AUC over 0.7 and 0.8, respectively. Average AUC was greater than 0.7 for all best-practices protein families demonstrating consistent enrichment performance across a broad range of proteins and ligand chemotypes. In both Astex and DUD datasets, docking performance is significantly improved employing a best-practices preparation scheme over using minimally-prepared structures from the PDB. Enrichment results for WScore, a new scoring function and sampling methodology integrating WaterMap and Glide, are presented for four DUD targets, hivrt, hsp90, cdk2, and fxa. WScore performance in early enrichment is consistently strong and all systems examined show AUC > 0.9 and superior early enrichment to DUD best-practices Glide SP results.

Water-assisted proton transfer in ferredoxin I.

The role of water molecules in assisting proton transfer (PT) is investigated for the proton-pumping protein ferredoxin I (FdI) from Azotobacter vinelandii. It was shown previously that individual water molecules can stabilize between Asp(15) and the buried [3Fe-4S](0) cluster and thus can potentially act as a proton relay in transferring H(+) from the protein to the µ(2) sulfur atom. Here, we generalize molecular mechanics with proton transfer to studying proton transfer reactions in the condensed phase. Both umbrella sampling simulations and electronic structure calculations suggest that the PT Asp(15)-COOH + H(2)O + [3Fe-4S](0) → Asp(15)-COO(-) + H(2)O + [3Fe-4S](0) H(+) is concerted, and no stable intermediate hydronium ion (H(3)O(+)) is expected. The free energy difference of 11.7 kcal/mol for the forward reaction is in good agreement with the experimental value (13.3 kcal/mol). For the reverse reaction (Asp(15)-COO(-) + H(2)O + [3Fe-4S](0)H(+) → Asp(15)-COOH + H(2)O + [3Fe-4S](0)), a larger barrier than for the forward reaction is correctly predicted, but it is quantitatively overestimated (23.1 kcal/mol from simulations versus 14.1 from experiment). Possible reasons for this discrepancy are discussed. Compared with the water-assisted process (ΔE ≈ 10 kcal/mol), water-unassisted proton transfer yields a considerably higher barrier of ΔE ≈ 35 kcal/mol.

Molecular Mechanics Force Field for Octahedral Organometallic Compounds with Inclusion of the Trans Influence.

Efficient calculation of the properties of metal-containing complexes relevant to catalysis is of major interest for better characterizing and optimizing the catalysts. For this, a new force field, called VALBOND-TRANS here, is proposed. It is based on the existing VALBOND force field of Landis and co-workers, extended by adding terms that account for electronic effects such as the trans influence of ligands on bond lengths and relative energies. Parameters and results for model octahedral complexes of Ru, Os, Rh, and Ir are determined and discussed. The model is then applied to the study of reactive intermediates involved in asymmetric hydrogenation catalyzed by iridium complexes with chiral phosphinooxazolines (PHOX) ligands. The new force field explores and capitalizes on the separation of electronic and steric effects on the stability of different diastereomers and reproduces DFT results which are consistent with experimental observations.

Search for non-nucleoside inhibitors of HIV-1 reverse transcriptase using chemical similarity, molecular docking, and MM-GB/SA scoring.

A virtual screening protocol has been applied to seek non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs) and its K103N mutant. First, a chemical similarity search on the Maybridge library was performed using known NNRTIs as reference structures. The top-ranked molecules obtained from this procedure plus 26 known NNRTIs were then docked into the binding sites of the wild-type reverse transcriptase (HIV-RT) and its K103N variant (K103N-RT) using Glide 3.5. The top-ranked 100 compounds from the docking for both proteins were post-scored with a procedure using molecular mechanics and continuum solvation (MM-GB/SA). The validity of the virtual screening protocol was supported by (i) testing of the MM-GB/SA procedure, (ii) agreement between predicted and crystallographic binding poses, (iii) recovery of known potent NNRTIs at the top of both rankings, and (iv) identification of top-scoring library compounds that are close in structure to recently reported NNRTI HTS hits. However, purchase and assaying of selected top-scoring compounds from the library failed to yield active anti-HIV agents. Nevertheless, the highest-ranked database compound, S10087, was pursued as containing a potentially viable core. Subsequent synthesis and assaying of S10087 analogues proposed by further computational analysis yielded anti-HIV agents with EC50 values as low as 310 nM. Thus, with the aid of computational tools, it was possible to evolve a false positive into a true active.

Elucidation of hydrolysis mechanisms for fatty acid amide hydrolase and its Lys142Ala variant via QM/MM simulations.

Fatty acid amide hydrolase (FAAH) is a serine hydrolase that degrades anandamide, an endocannabinoid, and oleamide, a sleep-inducing lipid, and has potential applications as a therapeutic target for neurological disorders. Remarkably, FAAH hydrolyzes amides and esters with similar rates; however, the normal preference for esters reemerges when Lys142 is mutated to alanine. To elucidate the hydrolysis mechanisms and the causes behind this variation of selectivity, mixed quantum and molecular mechanics (QM/MM) calculations were carried out to obtain free-energy profiles for alternative mechanisms for the enzymatic hydrolyses. The methodology features free-energy perturbation calculations in Monte Carlo simulations with PDDG/PM3 as the QM method. For wild-type FAAH, the results support a mechanism, which features proton transfer from Ser217 to Lys142, simultaneous proton transfer from Ser241 to Ser217, and attack of Ser241 on the substrate's carbonyl carbon to yield a tetrahedral intermediate, which subsequently undergoes elimination with simultaneous protonation of the leaving group by a Lys142-Ser217 proton shuttle. For the Lys142Ala mutant, a striking multistep sequence is proposed with simultaneous proton transfer from Ser241 to Ser217, attack of Ser241 on the carbonyl carbon of the substrate, and elimination of the leaving group and its protonation by Ser217. Support comes from the free-energy results, which well reproduce the observation that the Lys142Ala mutation in FAAH decreases the rate of hydrolysis for oleamide significantly more than for methyl oleate.

NO-MNDO: Reintroduction of the Overlap Matrix into MNDO.

The effect of reintroducing the overlap matrix into the secular equations for an NDDO (neglect of diatomic differential overlap)-based semiempirical molecular orbital method has been investigated. The modification is expected to improve the description of interactions between electron pairs. The idea has been tested by implementation and evaluation of a nonorthogonal version of the MNDO method (NO-MNDO) with parametrization for hydrogen, carbon, nitrogen, and oxygen. Overall, the accuracy of NO-MNDO for heats of formation is nearly identical to that for the more highly parametrized AM1 method. The mean absolute error (MAE) for heats of formation of a comprehensive set of 622 neutral, closed-shell molecules is reduced from 8.4 kcal/mol with MNDO to 6.8 kcal/mol with NO-MNDO. In addition, the performance for conformational equilibria and torsional barriers is significantly improved with NO-MNDO, presumably owing to the improved description of closed-shell interactions. For molecular geometries, the usual training and test sets have been expanded through use of MP2/6-31G(d) results for consistent comparisons. The performance of NO-MNDO for bond lengths, bond angles, and dihedral angles remains good with MAEs of 0.017 Å, 2.5°, and 4.5°. Additionally, NO-MNDO corrects severe errors by MNDO for R(•) + H-R' hydrogen-atom transfers, while testing for activation barriers for nine pericyclic reactions reveals only modest improvement.

Effects of Arg90 Neutralization on the Enzyme-Catalyzed Rearrangement of Chorismate to Prephenate.

Chorismate mutase (CM) is an enzyme that catalyzes the Claisen rearrangement of chorismate to prephenate. In a recent effort to understand the basis for catalysis by CM, Kienhöfer and co-workers (J. Am. Chem. Soc. 2003, 125, 3206-3207) reported results on the mutation of Arg90 in Bacillus subtilis CM (BsCM) to citrulline (Cit), an isosteric but neutral arginine analogue. An ca. 10(4)-fold decrease in kcat or 5.9 kcal/mol increase in the free-energy barrier (ΔG(‡)) for the overall catalysis was observed upon mutation. In this work, attention is turned to determining the key factors that contribute to the reduced catalytic efficiency of Arg90Cit BsCM. Using a combined QM/MM Monte Carlo/Free-Energy Perturbation method, a ΔΔG(‡) value of 3.3 kcal/mol is obtained. The higher free-energy barrier for the mutant is exclusively related to inferior stabilization of the TS, particularly one of its carboxylate groups, by neutral Cit. In addition, the reaction becomes 2.0 kcal/mol more exergonic. As BsCM is limited by product release, this step contributes to the remainder of the 10(4)-fold decrease in the rate constant in going from Arg90 to Cit.

Extension of the PDDG/PM3 Semiempirical Molecular Orbital Method to Sulfur, Silicon, and Phosphorus.

The PDDG/PM3 semiempirical molecular orbital method has been parameterized for molecules, ions, and complexes containing sulfur; the mean absolute error (MAE) for heats of formation, DeltaH(f), of 6.4 kcal/mol is 35 - 40 % smaller than for PM3, AM1, and MNDO/d. For completeness, parameterization was also carried out for silicon and phosphorous. For 144 silicon-containing molecules, the DeltaH(f) MAE for PDDG/PM3, PM3, and AM1 is 11 - 12 kcal/mol, while MNDO/d yields 9.4 kcal/mol. For the limited set of 43 phosphorus-containing molecules, MNDO/d also yields the best results followed by PDDG/PM3, AM1, and PM3. The benefits of the d-orbitals in MNDO/d for hypervalent compounds are apparent for silicon and phosphorous, while they are masked in the larger dataset for sulfur by large errors for branched compounds. Overall, for 1480 molecules, ions, and complexes containing the elements H, C, N, O, F, Si, P, S, Cl, Br, and I, the MAEs in kcal/mol for DeltaH(f) are 6.5 (PDDG/PM3), 8.7 (PM3), 10.3 (MNDO/d), 10.8 (AM1), and 19.8 (MNDO).

Extension of the PDDG/PM3 and PDDG/MNDO semiempirical molecular orbital methods to the halogens.

The new semiempirical methods, PDDG/PM3 and PDDG/MNDO, have been parameterized for halogens. For comparison, the original MNDO and PM3 were also reoptimized for the halogens using the same training set; these modified methods are referred to as MNDO' and PM3'. For 442 halogen-containing molecules, the smallest mean absolute error (MAE) in heats of formation is obtained with PDDG/PM3 (5.6 kcal/mol), followed by PM3' (6.1 kcal/mol), PDDG/MNDO (6.6 kcal/mol), PM3 (8.1 kcal/mol), MNDO' (8.5 kcal/mol), AM1 (11.1 kcal/mol), and MNDO (14.0 kcal/mol). For normal-valent halogen-containing molecules, the PDDG methods also provide improved heats of formation over MNDO/d. Hypervalent compounds were not included in the training set and improvements over the standard NDDO methods with sp basis sets were not obtained. For small haloalkanes, the PDDG methods yield more accurate heats of formation than are obtained from density functional theory (DFT) with the B3LYP and B3PW91 functionals using large basis sets. PDDG/PM3 and PM3' also give improved binding energies over the standard NDDO methods for complexes involving halide anions, and they are competitive with B3LYP/6-311++G(d,p) results including thermal corrections. Among the semiempirical methods studied, PDDG/PM3 also generates the best agreement with high-level ab initio G2 and CCSD(T) intrinsic activation energies for S(N)2 reactions involving methyl halides and halide anions. Finally, the MAEs in ionization potentials, dipole moments, and molecular geometries show that the parameter sets for the PDDG and reoptimized NDDO methods reduce the MAEs in heats of formation without compromising the other important QM observables.